* install scatterfit ado
net install scatterfit, from (https://raw.githubusercontent.com/leojahrens/scatterfit/master)

* install outreg2
ssc install outreg2

* open dataset (Indigenous_voices.dta)

* descriptive graphs
* Figure 1
sum pindi, det
scatterfit porce_rechazo pindi, fit(lfitci) opts(ylab(.42(0.1).95) xlab(0(5)90) name(fig1, replace)) coef 

* Figure 2
sum pcviolent, det
scatterfit porce_rechazo pcviolent, fit(lfitci) opts(ylab(.42(0.1).95) xlab(0(1)9.9) name(fig2, replace)) coef

* regressions
* with individual peoples
* OLS 
reg porce_rechazo c.pcviolent c.pmapuche c.prapanui c.paymara c.plican c.pquechua c.pcolla c.pdiaguita c.pkawesqar c.pyahgan c.pother c.CERCATOTAL c.kast2021 c.EDAD_MEAN c.lnpob c.pobreza c.educacion c.prural [aw=Cant_votos], vce(rob)
outreg2 using Appendix1.xls, replace  ctitle(OLS) alpha(0.001, 0.01, 0.05)
estimates store reg1
* GLM
glm porce_rechazo c.pcviolent c.pmapuche c.prapanui c.paymara c.plican c.pquechua c.pcolla c.pdiaguita c.pkawesqar c.pyahgan c.pother c.kast2021 c.CERCATOTAL c.EDAD_MEAN c.lnpob c.pobreza c.educacion c.prural [aw=Cant_votos], fam(bin) link(logit) vce(rob)
outreg2 using Appendix1.xls, append  ctitle(GLM) alpha(0.001, 0.01, 0.05)
estimates store reg2
* Figure 3
coefplot  (reg1, label(OLS)) (reg2, label(GLM)), keep(*pcviolent *pmapuche *rapanui *paymara *plican *pquechua *pcolla *pdiaguita *pkawesqar *pyahgan *pother) byopts(xrescale) xline(0) bycoefs name(fig3, replace)

* regressions with indigenous##rural interaction term
sum pindi, det
sum prural, det
reg porce_rechazo c.pindi##c.prural c.pobreza c.pcviolent c.CERCATOTAL c.kast2021 c.EDAD_MEAN c.educacion c.lnpob [aw=Cant_votos], vce(rob)
outreg2 using Appendix2.xls, replace alpha(0.001, 0.01, 0.05)
testparm c.pindi##c.pobreza
* margins at the range of values of indigenous population, and 1st, 5th, and 9th decile of rurality
margins, at(pindi=(2.1(5)87.8) prural=(10.81437 32.99587  72.42398 )) atmeans post
* Figure 4
marginsplot, name(fig4, replace) yscale(range(.4177449  .946957)) ylabel(.42(0.1).95, labsize(vsmall)) xlabel(0(5)100,labsize(tiny)) legend() title("") ytitle("Predicted Reject vote (%)") noci




